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We study the truncated Wigner method applied to a weakly interacting spinless Bose condensed 
£Nj ■ gas which is perturbed away from thermal equilibrium by a time-dependent external potential. 

The principle of the method is to generate an ensemble of classical fields ip(r) which samples the 
Wigner quasi-distribution function of the initial thermal equilibrium density operator of the gas, 
and then to evolve each classical field with the Gross-Pitaevskii equation. In the first part of the 
paper we improve the sampling technique over our previous work [Jour, of Mod. Opt. 47, 2629- 
' 2644 (2000)] and we test its accuracy against the exactly solvable model of the ideal Bose gas. 

In the second part of the paper we investigate the conditions of validity of the truncated Wigner 
method. For short evolution times it is known that the time-dependent Bogoliubov approximation is 
valid for almost pure condensates. The requirement that the truncated Wigner method reproduces 
the Bogoliubov prediction leads to the constraint that the number of field modes in the Wigner 

S simulation must be smaller than the number of particles in the gas. For longer evolution times the 
nonlinear dynamics of the noncondensed modes of the field plays an important role. To demonstrate 
this we analyse the case of a three dimensional spatially homogeneous Bose condensed gas and we 
test the ability of the truncated Wigner method to correctly reproduce the Beliaev-Landau damping 
^ ■ of an excitation of the condensate. We have identified the mechanism which limits the validity of 

the truncated Wigner method: the initial ensemble of classical fields, driven by the time-dependent 
Gross-Pitaevskii equation, thermalises to a classical field distribution at a temperature T c i ass which 
is larger than the initial temperature T of the quantum gas. When T c i ass significantly exceeds T a 
spurious damping is observed in the Wigner simulation. This leads to the second validity condition 
for the truncated Wigner method, T c i ass — T <C T, which requires that the maximum energy e ma x 
£^ ■ of the Bogoliubov modes in the simulation does not exceed a few ksT. 
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I. INTRODUCTION 



In Ref. [l[ the formalism of the Wigner representation of the density operator, widely used in quantum optics, 
was proposed as a possible way to study the time evolution of Bose-Einstein condensates in the truncated Wigner 
approximation Q- Like other existing approximate methods, such as the time-dependent Bogoliubov approach, it 

■ allows us to go beyond the commonly used Gross-Pitaevskii equation, in which the interactions between the condensate 
and the noncondensed atoms are neglected. Our aim in this paper is to illustrate the advantages and the limits of the 

■ truncated Wigner approach. 

For reasons of clarity we will address two different situations in two separate parts of the paper: (i) the case of a 
stationary Bose condensed gas in thermal equilibrium and (ii) a time-dependent case when the gas is brought out of 
equilibrium by a known external perturbation. Even if the stationary gas is the starting point for both situations, the 
problems raised by the application of the Wigner method are of a different nature in the two cases. 

(i) In the case of a Bose condensed gas in thermal equilibrium, the first step is to calculate the Wigner quasi- 
£j \ distribution function associated with the iV-body density operator a, which is a functional of a complex classical 
field ip(r). We showed in Q that this is possible in the Bogoliubov approximation when the noncondensed fraction 
of atoms is small. With such an approximation, the Hamiltonian of the system is quadratic in the noncondensed 
field and its Wigner functional is a Gaussian. After that, we went through some more technical work to calculate 
the Wigner functional of the whole matter field including the condensate mode. In our derivation we made further 
approximations in addition to the Bogoliubov approximation. This introduces some artifacts in the Wigner functional 
as far as the condensate mode is concerned p|. These artifacts are, however, insignificant when the number of thermally 
populated modes is much larger than one, or UbT S> Tuo in an isotropic trap of harmonic frequency u>, so that the 
fluctuations in the number of condensate particles, due to finite temperature, are much larger than one. Once the 
Wigner functional for the Bose condensed gas in thermal equilibrium is calculated, the goal is to be able to sample it 
numerically in order to compute averages of observables and probability distributions. In practice, this step consists 
in generating a set of random atomic fields {tp{r)} according to a probability distribution dictated by the Wigner 
functional. We have now developed a more efficient algorithm to sample the Wigner functional in the case of spatially 
inhomogeneous condensates in a trapping potential than the one that we had presented in a previous paper [|| , which 
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we will explain here in detail. As far as the equilibrium Bose condensed gases are concerned, our method in its regime 
of validity, is equivalent to the f/(l) symmetry-preserving Bogoliubov approach of @, @, up to second order in the 
small parameter of the theory, which is the square root of the noncondensed fraction. Compared with the traditional 
Bogoliubov approach, our method presents, however, the practical advantage of avoiding the direct diagonalisation 
of the Bogoliubov matrix, which is a heavy numerical task in 2D and 3D in the absence of rotational symmetry. 
Moreover, due to the stochastic formulation we adopt, our method gives us access to single realisations and to the 
probability distribution of some observables such as the number of condensate particles, not easily accessible by the 
traditional Bogoliubov method. We show some examples where we compare the probability distribution of the number 
of condensate particles obtained with our method with an exact calculation in case of the ideal Bose gas. 

(ii) Let us now consider the situation of a Bose condensed gas at thermal equilibrium which is brought out of 
equilibrium by a perturbation. The initial Wigner functional then evolves in time according to a kind of Fokker- 
Planck equation containing first and third order derivatives with respect to the atomic field. Numerical simulation 
of the exact evolution equation for the Wigner functional has intrinsic difficulties, as one would expect, since it 
represents the exact solution of the quantum many-body problem jjjj. We are less ambitious here, and we rely on an 
approximation that consists in neglecting the third order derivatives in the evolution equation. This is known as the 
truncated Wigner approximation [1]. For a delta interaction potential between a finite number of low energy modes 
of the atomic field, the third order derivatives are expected to give a contribution which is smaller than that of the 
first order derivatives when the occupation numbers of the modes are much larger than unity. If we reason in terms of 
the stochastic fields ip(r, t) which sample the Wigner distribution at time t, then the truncated Wigner approximation 
corresponds to evolving the initial set of stochastic fields according to the Gross-Pitaevskii equation Q: 



.^A + U(r,t)+g\i,\ 2 
2m 



V>, (i) 



where r is the set of single particle spatial coordinates, m is the atom mass, U is the trapping potential and g is 
the coupling constant originating from the effective low energy interaction potential V(r\ — r%) — gS(r\ — r%) and 
proportional to the s-wave scattering length a of the true interaction potential, g — 47r?i 2 a/m. Here, the crucial 
difference with respect to the usual Gross-Pitaevskii equation is that the field is now the whole matter field rather 
than the condensate field. 

This procedure of evolving a set of random fields with the Gross-Pitaevskii equation is known as the classical field 
approximation, since equation ([1]) can be formally obtained from the Heisenberg equation of motion for the atomic 
field operator tp by replacing the field operator by a c-number field. The classical field approximation has already 
been used in the Glauber-P representation to study the formation of the condensate [jj EH EH, E2, EH ■ We face here a 
different situation: we assume an initially existing condensate and we use the Wigner representation, rather than the 
Glauber-P representation. The Wigner representation is in fact known in quantum optics to make the classical field 
approximation more accurate than in the Glauber-P representation because the "right amount" of quantum noise is 
contained in the initial state [Llj . For a single mode system with a Kerr type nonlinearity and an occupation number 
n, the term neglected in the Wigner evolution equation is a third order derivative which is 1/n 2 times smaller than 
the classical field term, whereas the term neglected in the Glauber-P evolution equation is a second order derivative, 
which is only 1/n times smaller than the classical field term. In the case of Bose-Einstein condensates however, we 
face a highly multimode problem and, therefore, the accuracy of the truncated Wigner approach needs to be revisited. 
We approach this problem in the second part of the paper. The strategy we adopt is to compare the predictions of the 
truncated Wigner method with existing well-established results: first with the time-dependent Bogoliubov approach 
and then with the Landau-Beliaev damping of a collective excitation in a spatially homogeneous condensate. 



II. BASIC NOTATIONS AND ASSUMPTIONS 
A. Model Hamiltonian on a discrete grid 

Let us express a simple quantity like the mean atomic density using the Wigner representation: 

(^(r)^(r)) = (p(r)1>(r)) w - ~{$(r),ft(r)]), (2) 

where (. ..}w represents the average over the Wigner quasi-distribution function. This shows that the discretisation 
of the problem on a finite grid is necessary to avoid infinities: in the continuous version of the problem, [V>(r), ^'{f)) = 
6(0) = +oo. Physically this divergence comes from the fact that, in the Wigner point of view, some noise is included 
in each mode of the classical field ip to mimic quantum noise; this extra noise adds up to infinity for a system with 



3 



an infinite number of modes. Therefore we use, from the beginning, a discrete formulation of our problem which will 
make it also suitable for numerical simulations. 

We consider a discrete spatial grid forming a box of length L v along the direction v — x,y, z with an even number 
n v of equally spaced points. We denote M = \\ u n v the number of points on the grid, V = Y[ u L, y the volume of the 
grid and dV = VjM the volume of the unit cell of the grid. We take periodic boundary conditions in the box [ill ]. 
We can then expand the field operator over plane waves 

V-(r)-^a fc 4r e * r ' (3) 



where au annihilates a particle of momentum k and where the components of k are fc„ = 2itj v / L v with the integers 
j v running from — n„/2 to n v /2 — 1. We then have the inverse formula: 

a k = dVj2^- ihr i>(r)- (4) 
For each node r% on the spatial grid, we find the commutation relations for the field operator: 

mn),4>Hr j )] = -L5 i , j (5) 
and the discretised model Hamiltonian that we use is: 

H = ^~4 a * + dV J2 U(r)4>Hr)4>(r) + |dF^^(r)^( r )^( r )^( r ) . ( 6 ) 

k r r 

The first term in ([5]) is the kinetic energy, which is easy to calculate in the momentum representation. In the position 
representation, the kinetic energy is a matrix that couples the Af points of the grid. In the following we will write it 
as p 2 /2m. The second term is the trapping potential. The last term represents the atomic interactions modeled by a 
discrete Kronecker 5 potential 

V(ri -r 2 ) = — 5 ri ,r 2 , (7) 

with a coupling constant g — 4:nh 2 a/m, where a is the s-wave scattering length of the true interaction potential. 

We indicate briefly some requirements for the discrete Hamiltonian to be a good representation of reality. First, 
the spatial step of the grid should be smaller than the macroscopic physical scales of the problem: 

dx v -C £ and dx v <C A, (8) 



where £ = l/y^87rp|af is the healing length for the maximal atomic density p and A = \l 2tt?i 2 /mksT is the thermal 
de Broglie wavelength at temperature T. Secondly, the spatial step of the grid should be larger than the absolute 
value of the scattering length a: 

dx v ^>\a\, (9) 

so that the scattering amplitude of the model potential is indeed very close to a. Another way of saying this is 
that the model potential {7} can be treated in the Born approximation for the low energy waves. A more precise 
treatment, detailed in the appendix [Aj is to replace in ([7|) the coupling constant g by its bare value go adjusted so 
that the scattering length of the model potential on the grid is exactly equal to a. 



B. Wigner representation 



The Wigner quasi-distribution function associated with the TV-body density operator a is defined as the Fourier 
transform of the characteristic function x' 



W(ip) 



J -q dRe~f(r)dlm-y(r)dV e dv Y, r 7"»V>M-7W(r) 



X(7) = Tr \ae 



dV ^2 ^{r)^ {r)— 7* (r)i/>(r) 



(10) 

(ii) 
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where j(r) is a complex field on the spatial grid and a is the density operator of the system. With this definition the 
Wigner function is normalised to unity: 



Y[dRei>(r)dIm^(r)dV W(ip) = 1. 



(121 



We recall that the moments of the Wigner function correspond to totally symmetrised quantum expectation values, 
i.e., 



{Ox 



, O n ) W = ^2 Tr [Op(l) ■ ■ ■ Op(r 



(13) 



where the sum is taken over all the permutations P of n objects, Ok stands for ip or ip* in some point of the grid and 
Ofc is the corresponding field operator. 

The equation of motion for the density operator a 



d 1 r - . 

dt a = m [H > a] 

can be written exactly as the following equation of motion for the Wigner distribution: 



in—— = > 
dt ^ 







with a drift term 



dip(r) 



(-UW) 



a 3 



A{dV) 2 d 2 ip{r)dip*(r) 



(ip{r)W) - c.c. 



f- + tf(r,t) + <^V- 4j 
2m dv 



(14) 



(15) 



(16) 



The truncated Wigner approximation consists in neglecting the cubic derivatives in the equation for W. The resulting 
equation reduces to the drift term whose effect amounts to evolving the field ip according to an equation which 
resembles the Gross-Pitaevskii equation (|T|). The constant term —g/dV inside the brackets of the above equation can 
be eliminated by a redefinition of the global phase of tj), which has no physical consequence for observables conserving 
the number of particles. 



III. SAMPLING THE WIGNER FUNCTIONAL FOR A BOSE CONDENSED GAS IN THERMAL 

EQUILIBRIUM 

In [|[ we derive an expression of the Wigner functional for a Bose condensed gas in thermal equilibrium in the 
frame of the U(l) symmetry-preserving Bogoliubov approach 0, @, in which the gas has a fixed total number of 
particles equal to N. We first introduce the approximate condensate wavefunction 4>{r), which is a solution of the 
time-independent Gross-Pitaevskii equation: 



H gp cj) = 



2-+U(r,t = 0) + N 9 \<t>\ 1 
2m 



I* 



= 0. 



(17) 



We then split the classical field tjj(r) into components orthogonal and parallel to the condensate wavefunction <fi(r): 

i/f(r) = a^(r)+yj ± (r) (18) 
a = dVY,<f>*(rMr). (19) 

r 

The Wigner functional provides us with the joint probability distributions of the transverse classical field ijjj_(r), that 
we call the noncondensed field, and of the complex amplitude a^. Due to the U(l) symmetry-preserving character of 
the theory, the final Wigner functional is of the form [3j 



dd 
2^ 



Wo(e- iS V0- 



(20) 



This means that one can sample the distribution W(tf>) by (i) choosing a random field t/j according to the distribution 
Wo(tp)> (ii) choosing a random global phase 9 uniformly distributed between and 27r, and (iii) forming the total 
atomic field as ^tot( r ) = e l9 -0(r). In practice, the global phase factor e %e is unimportant to calculate the expectation 
value of observables that conserve the number of particles. Since the other observables have a vanishing mean value, 
we can limit ourselves to the sampling of the 9 = component of the Wigner functional, Wo(ip). 
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A. Sampling the distribution of the noncondensed held 

The first step of the sampling procedure consists in generating a set of noncondensed fields {ip±} according to the 
probability distribution 



P(ipx) ot exp 



-dV WlM-M^ 



(21) 



where we have collected all the components of ipx and ip^_ in a single vector with 2TV components, M is the 27V x 27V 
matrix: 



M = r) tanh 



2k B T 



with 



V = 



1 

-1 I ' 



and where £ is a 27V x 27V matrix, which is the discretised version of the Bogoliubov operator of || : 



£ = 



H ep + NgQ\0\ 2 Q NgQ(f> 2 Q* 
-NgQ*4>* 2 Q -H* p -NgQ*\4>\ 2 Q* 



(22) 



(23) 



(24) 



In this expression the TV x TV matrix Q projects orthogonally to the condensate wavefunction <p in the discrete spatial 
grid {r{\ representation, 



Qij = Sij - dVct>(ri)4>*(rj). 
Note that the matrix M can be shown to be Hermitian from the fact that & = rjCij. 



(25) 



1. Direct diagonalisation of C 



If the eigenvectors of C are known, we can use the following modal expansion: 



V>1 



5> 



K ( v i 



(26) 



where the sum is to be taken over all cigenmodcs (uk,Vk) of £ normalisable as (Uk\Uk) — (wfc|ufc) = 1, with 
corresponding eigenvalues Since the condensate is assumed to be in a thermodynamically stable or metastablc 
state, all the ek are positive [1 61 ] . The probability distribution (|21[) is then a simple product of Gaussian distributions 
for the complex amplitudes bk- 



Pk(b k ) 



' tanh ( 'oT^ 7 f 
7T \2kb-1 



exp 



-2|6 fc | 2 tanh 



\2k B T 



(27) 



Each Gaussian distribution is easily sampled numerically [17| . Note that, even at zero temperature, the Gaussian 
distribution has a nonzero width: this is a signature of the extra noise introduced in the Wigner representation to 
mimic quantum noise. 



2. Brownian motion simulation 



The sampling of the distribution (|21| can actually be performed without diagonalisation of C (an advantage for 
spatially inhomogeneous Bose condensates in the absence of rotational symmetry [4|) by means of a Brownian motion 
simulation for the noncondensed field: 



d\ t)=-adt(^ 



Y 



di 

da* 



(28) 
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where the field e?£ is the noise term. The time t here is a purely fictitious time with no physical meaning and will be 
taken to be dimensionless. On our discrete grid, ?p± is a vector with TV components, d£ is a Gaussian random vector 
of TV components with zero mean and a covariance matrix (d&dt;*) equal to (2dt / dV)5ij , while a, Y are 27V" x 27V 
matrices. To ensure that the Brownian motion relaxes towards the correct probability distribution (I21[) we require 
that the drift matrix a and the diffusion matrix D = Y(Y<) satisfy a generalised Einstein's relation [4|: 

D~ 1 a = a i D- 1 = 2M, (29) 

where M is the matrix (|22| . There is, of course, no unique choice for a and Y. With respect to our previous work 
we have largely improved the efficiency of our simulation by a different choice of a, Y and by the use of a second 
order integration scheme of the stochastic differential equation (|28p. more efficient than the usual first order Euler's 
scheme. In the appendix[B]we give a detailed description of these improvements, useful to the reader who is interested 
in implementing the numerical algorithm. 



B. Sampling the condensate amplitude 

We now have to sample the condensate amplitude from the Wigner functional Wq . This amplitude turns out to 
be real, and can be written as 

= v^Vo where N = a^aj, . (30) 

Since we already know how to generate the noncondensed part of the field -ipj_, we have to sample the conditional 
distribution P(N \ip±). 

Due to a first approximation that we have performed in Q , which consists in treating "classically" the condensate 
mode and neglecting its quantum fluctuations in the limit of a large number of condensate particles, the probability 
distribution P(Nq), that we will obtain by averaging P(No\ip±) over the stochastic realisations of the noncondensed 
field ip±, actually coincides with the probability distribution of the number of condensed particles aXa^ so that within 
this approximation we have: 

(N ) = (fifa), (31) 

Var(iV ) = Var(a^),... (32) 

Note that this should not be the case for the exact Wigner distribution as, e.g., the average (Nq) should be equal to 
(a^a^) + 1/2 and the variance of No should exceed the variance of a^a,/, by 1/4. 

We show in 3] that, when the number of thermally populated modes is much larger than one, the width in No of 
the conditional distribution P(No\ipx) is much narrower than the width of the distribution P(Nq), so that we can 
replace the distribution P(No\ip±) by a delta function centered on its mean value. With this second, more severe, 
approximation we get for the sampling: 

7V ^ Mean(iV |^) = C - \dV^lM ' [H - M 2 ] (t*) , (33) 

where the constant C is finite only in the discretised version and is given by 

C = N - -Tr M + -Tr Q. (34) 
4 2 

Here, the trace of the projector Q is simply the number of modes in the simulation minus one. 

The second approximation (J33J) does not introduce errors in the average (No). We are able to verify a posteriori that 
the error introduced in the variance (Nq) — (Nq) 2 is small in the following way: on one hand we calculate the variance 
of iVo (Var(iVo)), by using (|33|) . On the other hand we calculate the variance V&r(5N) of the number of noncondensed 
particles by using directly the ensemble of noncondensed fields {f/'J-}- Since the total number of particles is fixed one 
should have Var(7V ) = Var^a^) = Var^^-^jJ, and deviation from this identity gives us the error of Var(7V ). 

We are now ready to form the total field: 

f cb (2) (r)\ 

i>(r) = ^oU(r) + ^ 1 \+ip±(r). (35) 

(2) 

The function (j> ± is a correction to the condensate wavefunction including the condensate depletion neglected in 
the Gross-Pitaevskii equation (|17p and the mean field effect of the noncondensed particles. This correction can be 
calculated from the ensemble of noncondensed fields a s explained in [J]. As we will see in section [IV Al its 

contribution to the one-body density matrix is of the same order as that of ip± and therefore has to be included. 
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C. Tests and applications: Distribution of the number of condensate particles 



We can use the sampling procedure described above to calculate some equilibrium properties of the Bose condensed 
gas. Recently, the variance of the number of particles in the condensate has drawn increasing attention [l8l [T9L [20|. 
In our case we have access to the whole probability distribution for N by applying equation (|33p to the ensemble of 
stochastic noncondensed fields {il>±}- 



1. Ideal Bose gas 



As a test we check our probability distribution for the number of condensate particles against the exact one for the 
ideal Bose gas (g — 0) in one and two dimensions. The results are in figure [TJ 
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FIG. 1: Probability distribution in the canonical ensemble of the number of condensate particles for the ideal Bose gas in 
thermal equilibrium in an isotropic harmonic potential U(r) = \rmJ 2 'r 2 . (a) In a ID model for feT = 30huj, and N = 10 000. 
For the Wigner simulation 2000 realisations have been performed on a grid with 128 points. For the exact Bogoliubov rejection 
method described in the end of this subsection on the ideal gas, 400 000 realisations have been performed so that the statistical 
error is less than one per cent for the most populated channels of the histogram, (b) In a 2D model for ksT — 307kj, and 
N = 8 000. For the Wigner simulation 500 realisations have been performed on a grid with 128 x 128 points. For the exact 
sampling 100 000 realisations have been performed. 



The distributions of the number of condensed particles iVo are clearly not Gaussian. To characterise them, besides 
the mean and the variance of Nq one can introduce the skewness defined as: 

skewW = M^m- (36) 

For the parameters of figure [1] we give the mean, the standard deviation and the skewness of Nq obtained from the 
simulation, together with their exact values: 





ID simulation 


ID exact 


2D simulation 


2D exact 


(No) 


9882. 


9880. 


6403. 


6415. 


ATVo 


37.5 


38.3 


75.9 


77.1 


skew(iVo) 


-1.20 


-1.16 


-0.40 


-0.334 



In what follows we explain in some detail how the exact probability distribution for the ideal Bose gas is obtained. 
Let a be the density operator for the ideal Bose gas in the canonical ensemble: 

&=^e- f>A p N . (37) 



8 



The operator p^ projects onto the subspace with N particles, and H = ^ fc e k a\a k is written in the eigenbasis of the 
trapping potential. In the spirit of the number conserving Bogoliubov method, we eliminate the condensate mode by 
writing 

a\ao = N - a\a k . (38) 

Since the total number of particles is fixed we can replace the operator N by the c- number N in (|38[) . Furthermore we 
establish a one to one correspondence between (i) each configuration of excited modes {n k , k > 0} having a number 
of excited particles N' = n k lower than TV and (ii) each configuration of the whole system with particles in 
excited mode k and N — N particles in the condensate. We then obviously have to reject the configurations of excited 
modes for which the number of particles in the excited states TV' is larger than N. This amounts to reformulating the 
effect of the projector p^ in terms of an Heaviside function Y . We then rewrite a as: 

^ = I e -/3eoJV e -/5E^- e °) a l a *y In-^&I&u ) . (39) 



z 



For the sampling procedure we use a rejection method i.e. we sample the probability distribution of the number of 
particles n k in each mode k ^ without the constraint imposed by the Heaviside function and we reject configurations 
with TV' > N. In this scheme we have to generate the rife, k = 1, . . . , Af, according to the probability distribution 

Pk(n k )=\Z k (l-\ k ) with Afe=e-^- £ °). (40) 

For each k we proceed as follows: in a loop over n k starting from we generate a random number e uniformly 
distributed in the interval [0,1] and we compare it with X k : if e < Afc, we proceed with the next step of the loop, 
otherwise we exit from the loop and the current value of n k is returned. 

The calculation can also be done in the Bogoliubov approximation, that is by neglecting the Heaviside function 
in (f39|) . For the parameters of figure [T] this is actually an excellent approximation, as the mean population of the 
condensate mode is much larger than its standard deviation, and the corresponding approximate results are in practice 
indistinguishable from the exact ones. The predictions of this Bogoliubov approximation for the first three moments 
of iVo involve a sum over all the excited modes of the trapping potential: 



(No) = 



fc#0 



Var(TV ) = ^n fe (l + n k ) 



k^O 

((N Q - (N )f) = J2 2n l + 3n| + n k (41) 

k^O 

where n k = l/(exp(/3(efe — eo)) — 1) is the mean occupation number of the mode k. In the limit fc^T 3> hu> for an 
isotropic harmonic trap an analytical calculation, detailed in the appendix [Cj shows that the skewness tends to a 
constant in ID, tends to zero logarithmically in 2D and tends to zero polynomially in 3D plj : 

skew 1D (TV ) ~ -^S§2 =-1-139547... 



skew 2 D(-/V ) 



2(C(2)+C(3)) 



(\og(k B T/hcj) + 1 + 7 + C(2)) 3 / 2 



k (N) log(WM+7+| + 3C(2) + 2C(3) 

s ew 3Di oj (fc B T/M 3/2 {C(2) + (3hLu/2k B T)[log(k B T/hu;) + 7 + 1 - C(2)/3]} 3 / 2 1 ' 

where £ is the Riemann Zeta function and 7 = 0.57721 ... is Euler's constant. 



2. Interacting case 



As an example we show in figure [2] the probability distribution for the number of condensate particles in the 
interacting case to demonstrate that the large skewness of No in ID can even be enhanced in presence of interaction: 
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the skewness of Ao in figure [2] is equal to —2.3. We have been able [22| to calculate P(Nq) in the Bogoliubov 
approximation in the interacting case starting from the sampling of the Wigner distribution of the noncondcnscd field 
(|2"Tj) . We compare the results with the Wigner approach in the same figure. As expected the agreement is excellent 
in the regime UbT = 30fko 3> %u>. 



Bogoliubov 

Wigner 





















8000 8500 9000 9500 10000 



FIG. 2: Probability distribution of the number of condensate particles in the canonical ensemble for a ID interacting Bose gas 
in thermal equilibrium in a harmonic trap U{x) = imw 2 i 2 , with fcsT = SOhcu, /i = U.lhu) and N = 10 000, corresponding to 
a coupling constant g = 0.0lhu}(h/mu>) 1 ^ 2 . The results have been obtained with the Wigner method using 2000 realisations on 
a grid with 128 points. The dashed line is the histogram of the probability distribution of No in the Bogoliubov approximation 
generated using the same 2000 realisations, obtained with a method described in 



IV. THE TRUNCATED WIGNER METHOD FOR A TIME-DEPENDENT BOSE CONDENSED GAS 



In this section we investigate the conditions of validity of the truncated Wigner approach for time-dependent Bose- 
Einstcin condensates. The strategy that we adopt is to compare the predictions of the truncated Wigner approach to 
well-established theories: the time-dependent Bogoliubov approach in section IIV Al and the Landau-Beliaev damping 
of a collective excitation in a spatially homogeneous condensate, in section HVBI 



A. The truncated Wigner method vs the time-dependent Bogoliubov method 

In this section we investigate analytically the equivalence between the time-dependent Bogoliubov approach of Q 
and the truncated Wigner method in the limit in which the noncondensed fraction is small. 

We begin by sketching the number conserving Bogoliubov method of Ref. || . We split the atomic field operator 
into components parallel and orthogonal to the exact time-dependent condensate wavefunction (j> cx [23[ (omitting for 
simplicity the time label for the field operators and for the condensate wavefunction): 

$ ( r ) = «0c X 0c X (r) + ip± (r) (43) 

and we consider the limit 

N — > oo N g — constant T = constant Af = constant. (44) 

In 1 one performs a formal systematic expansion in powers of 1/ y~N of the exact condensate wavefunction cf> cx 

4>Ur) = t(r) + + t-LL + . . . (45) 



and of the noncondensed field 



A. 



cx (r) = ^=al V±(r) = A(r) + ^=A« (r) + . . . . (46) 



AT * e * VN 
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Note that in the lowest order approximation to A ex the exact condensate wavefunction CX is replaced by the solution 
<j) of the time-dependent Gross-Pitaevskii equation 

ihd t (j) = [p 2 /2m + U(r, t) + Ng\<j)\ 2 } <j> (47) 
and a^/v^/V is replaced by the phase operator — a^aXa^) so that 



A(r) = -^=4 ^) ~ 4>{r)h (48) 



and A(r) satisfies bosonic commutation relations 



[A(r),At( s )] = — Q r , s (49) 



1 

dV' 

where the matrix Q r . s = 5 r , s — dV(b{r)(b* (s) projects orthogonally to (b. To the first two leading orders in 1/yN one 
obtains an approximate form of the one-body density matrix: 

(r\p\a) = (ft(8)$(r)) = (N-{5N))cb(r)ct>*(s) 

+ (At(a)A(r)) 

+ r(s)0f(r)+0(r)0f*(s) 

+ (50) 



We call the first term "parallel-parallel" because it originates from the product of two parts of the field both parallel 
to the condensate wavefunction; it describes the physics of a pure condensate with N — (8N) particles. The second 
term, which we call "orthogonal-orthogonal" because A is orthogonal to <fi, describes the noncondensed particles in the 
Bogoliubov approximation. The third term, called "orthogonal-parallel" , describes corrections to the Gross-Pitaevskii 
condensate wavefunction due to the presence of noncondensed particles [5|. In (|50p (SN) is the average number of 
noncondensed particles in the Bogoliubov approximation: 

(SN) =J2 dV(AHr)A(r)}. (51) 

r 

The evolution equations for A and (b^ are given in appendix |Dj 

Having described the Bogoliubov method, let us now consider the truncated Wigner approach in the limit (|4"4"|) . We 
expand the classical field in powers of 1 / vN : 

ib = VNib^ + + -Uv (2) + • ■ ■ (52) 



where the tpV' are of the order of unity. We immediately note that the leading term of this expansion corresponds to a 
pure condensate with N particles so that ib^ is simply the solution of the time-dependent Gross-Pitaevskii equation 
P?)) . V'*' '* = <t>- This physically clear fact will be checked explicitly in what follows. In the initial thermal equilibrium 
state at time t = we expand (pj|) in powers of 1/VN : 



so that we can identify explicitly: 



= VN - SN = VN - -^L + ... (53) 



2 JN 



^°)(i 


= 0) = 






= 0) = v± 




^ 2 \t 







(54) 
(55) 

(56) 
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Following the same procedure as in the quantum case, we split each term of the expansion into a component along <p 
and a component orthogonal to (\>: 

^U)(r) =^')^(r)+V$' ) . (57) 

We calculate now the one-body density matrix p. Since we are using the Wigner representation for the atomic field 
on a finite spatial grid we have: 

(r\p\s) = <V»tf(r)> - ^ r , s (58) 

where dV is the unit cell volume of the spatial grid and S riS is a Kronecker S. Note that to simplify the notation we 
have omitted the subscript W on the right hand side of the equation since the quantum and Wigner averages can be 
readily distinguished by the hats on the operators. We insert the expansions ([52"]) and ([57)) into ([55)) and we use the 
fact that ip^ — 4> to obtain: 



{r\p\s) T w = (j>*{s)(j>{r) 



n + Viv<e (1) + + (ie (1) i 2 ) + <e (2) + e (2) *) - 1 



+ ^>{ s )^l\r))-^Q r , s 

+ ct>*{s)[^N{^lHr)) + (e (1) >i 1} (r)) + (Vf(r))] + {r «- s}* 

+ ( 59 ) 



N 

where we have collected the terms "parallel-parallel" in the first line, the terms "orthogonal-orthogonal" in the second 
line and the terms "orthogonal-parallel" in the third line, and where the matrix Q r , s /dV = & r . s /dV — <fi(r)(f>* (s) is the 
discrete version of the projector Q = 1 — |</>)(</>|. As we show in appendix[El by using the evolution equation of the 
field ([T]) and the initial conditions ([5~4")) , (|55p and (1551) the following identities hold at all times: 

</> (0) = 4> (60) 

iv(e (1) +^ (1) *) + (ie (1) i 2 ) + (e (2) +^ (2) *) = -(6N) (6i) 

(^*(s)^l\r))-^Q r , s = (At(«)A(r)) (62) 
N(^l\r)) + (^*^\r)) + (^\r)) = ^(r). (63) 



As we have already mentioned the first identity l|60p reflects the fact that at zero order in the expansion we have a 
pure condensate with TV particles evolving according to the time-dependent Gross-Pitaevskii equation. At time t = 
the three other identities are easily established since we have simply (V^) = 0, = and £ (2) = -5N/2. At 

later times the mean value remains equal to zero while £^ develops a nonzero imaginary part corresponding 

to phase change of ?/> in the mode <fi due to the interaction with the noncondcnscd particles 

tp = \/N<f> + + . . . ~ VNe^ 1)/VW <j> +... (64) 

After averaging over all stochastic realisations, this random phase change contributes to the condensate depletion in 
(IBT)) and to the correction (f>^ to the condensate wavefunction in ([6"3")l [24j ] . As a consequence of the purely imaginary 
character of the quantity proportional to y/~N in (|6ip vanishes. The identity ([B"2")) reflects the fact that in the 
linearised regime quantum fluctuations (here A) and classical fluctuations (here ipi^) around the Gross-Pitaevskii 

condensate field sfNcj), evolve according to the same equations. We find interestingly that the average {^±^) m 
([63)) evolves under the influence of the mean field of the noncondensed particles, i.e. the Hartree-Fock term and the 



anomalous average contribution. In the Wigner representation the Hartree-Fock mean field term 2g(ip^*ip ( ^') differs 
from the physical mean field 2g(A j 'A} by the term ^(1 — \</)\ 2 dV)/dV ~ g/dV. We note however that this brings in 
a global phase change of the condensate wavefunction having no effect on the one-body density matrix, and which is 
compensated anyway by the —g/dV term in the Wigner drift term (|16|) . In our calculations this is reflected by the 
fact that this term does not contribute to <p]_ ■ 
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With the identities (|60H63p we identify line by line the quantum expression (|50| and the truncated Wigner expression 
(f59"| for the one-body density matrix of the system up to terms of 0(1): these two expressions coincide apart from 
the term 1/2 in the occupation number of the mode <fi. This slight difference (1/2 <C N) comes from the fact that in 
the initial sampling of the Wigner function in thermal equilibrium we have treated classically the condensate mode. 
These results establish the equivalence between the truncated Wigner method and the time-dependent Bogoliubov 
approach of || up to neglected terms 0(1/ y/N) in the limit (J3H). 

Let us however come back to the expansions performed in the limit (I44|) . We have mentioned that the small formal 
parameter is 1/y/N but we now wish to identify the small physical parameter of the expansion. In the quantum 
theory of [5J] one gets the small parameter 



Equant 



- \ 1/2 

(5N) \ 



(65) 



where (SN) is the Bogoliubov prediction for the number of noncondensed particles. In the expansion (|52[) of the 
evolving classical field we compare the norm of the first two terms, ignoring the field phase change 



'{dV^ r \^\ 2 )\ 1/2 _ f(5N) + (Af-l)/2^ ' 2 



[■■Wig 



(66) 



N J V N I 

The validity condition of the expansion (|52[) in the truncated Wigner approach is then: 

N > (SN) , Af/2 (67) 

which is more restrictive than in the quantum case. What indeed happens in the regime (SN) <C N < Af/2? We 
expect the truncated Wigner approach not to recover the predictions of the Bogoliubov approach of [f| which are 
correct in this limit. We therefore set a necessary condition for the validity of the truncated Wigner approach: 

iV»jV/2. (68) 

We interpret this condition as follows: the extra noise introduced in the Wigner representation (see discussion after 
([2"7])) contributes to the nonlinear term g\ifj\ 2 in the evolution equation for the field; (|6"8|) means that this fluctuating 
additional mean field potential of order g/(2dV) should be much smaller than the condensate mean field of order 
gN/V where V = AfdV is the volume of the system. Condition (f6"B")) is also equivalent to pdV ^> 1, where p is the 
atomic density, i.e. there should be on average more than one particle per grid site. We note that it is compatible with 
the conditions ([5]) on the spatial steps of the grid in the regime of a degenerate (pA 3 3> 1) and a weakly interacting 
(p£ 3 3> 1) Bose gas. Condition ([68]) is therefore generically not restrictive. 

A last important point for this subsection is that the time-dependent Bogoliubov approach, relying on a linearisation 
of the field equations around a pure condensate solution, is usually restricted to short times in the case of an excited 
condensate, so it cannot be used to test the condition of validity of the truncated Wigner approach in the long time 
limit. It was found indeed in [25[ that nonlinearity effects in the condensate motion can lead to a polynomial or 
even exponential increase in time of (SN) which eventually invalidates the time-dependent Bogoliubov approach. The 
truncated Wigner approach in its full nonlinear version does not have this limitation however, as we have checked 
with a time-dependent ID model in Q. 

B. Beliaev-Landau damping in the truncated Wigner approach 

In this section we consider a spatially homogeneous Bose condensed gas in a cubic box in three dimensions with 
periodic boundary conditions. We imagine that with a Bragg scattering technique we excite coherently a Bogoliubov 
mode of the stationary Bose gas, as was done experimentally at MIT [26J, |27J , and we study how the excitation decays 
in the Wigner approach due to Landau and Beliaev damping. 

1. Excitation procedure and numerical results 

We wish to excite coherently the Bogoliubov mode of wavevector ko ^ 0. With a Bragg scattering technique using 
two laser beams with wave vector difference q and frequency difference u> we induce a perturbation potential 

W= f d 3 r (^ e z( - q r -^ +c.c) (69) 
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We match the wavevector and frequency of the perturbation to the wavevector kg and the eigenfrequency loq — sq/K 
of the Bogoliubov mode we wish to excite: 

q = fc oj = e /h = lj . (70) 

During the excitation phase, we expect that two Bogoliubov modes are excited from the condensate, the modes with 
wavevectors kg and —kg. We anticipate the perturbative approach of next subsection which predicts that the mode of 
wavevector ko, being excited resonantly, has an amplitude growing linearly with time, while the mode with wavevector 
— fcoj being excited off-resonance, has an oscillating amplitude vanishing periodically when t is a multiple integer of 
tt/luq. In the truncated Wigner simulation we therefore stop the excitation phase at 

*exc = — • (71) 

We introduce the amplitudes of the classical field ip of the Bogoliubov modes. We first define the field 

A s tatic(r) = -=ali/j±(r) (72) 
V iV 



where and ipj_ are the components of ip orthogonal and parallel to the static condensate wavefunction <f)(r) — 1/L 3 / 2 
(see (fl8|) ). The component along the Bogoliubov mode with wavevector k is then 

b k = dV ]T u* k (r) A statlc (r) - v* k (r) A* atlc (r) . (73) 

r 

The functions itk and Vk are plane waves with wavevector k ^ 

u k (r) = --L U k e lk - r v k (r) = ^V k e lkr (74) 



1 / h 2 k 2 /2m \ 1/4 
U k + V k = — = J — (75) 



and the real coefficients U k and V k are normalised to U% — V k 2 = 1 

U k -V k ~ ijL A k 2 /2m + 2n 

where the chemical potential is fi = gN/L 3 . 

We denote by bo the amplitude of the field A sta tic along the Bogoliubov mode of wavevector ko, and 6_o the 
amplitude along the mode with opposite wavevector. We show the mean values of these amplitudes as function of 
time obtained from the truncated Wigner simulation in figure [3] In the initial thermal state these mean values vanish, 
and they become nonzero during the excitation phase due to the coherent excitation procedure. At later times they 
decay to zero again {28] . 



2. Perturbative analysis of the truncated Wigner approach: Beliaev-Landau damping 

In the appendix [F] we report the exact equations of motion of the classical field A sta ti c defined by (f72")) in the 
truncated Wigner approach. We now make the assumption that A sta tic is small compared with y/~N<f>, implying that 

N » (SN) , - (76) 

where (SN) represents here the mean number of particles in the excited modes of the cubic box. In this regime we 
neglect terms which are at least cubic in A sta tic in (|F2[) and we replace the number of particles in the ground state 
of the box by the total number of particles N, except in the zeroth order term in A stat j c where we replace it by its 
initial mean value (iVo). We then find: 

ih^_A static ~ y/(N ) Qh <jy + Q/i A static + (A* tatic + 2A static ) 

— ■ Q( A s t a tic A s t a tic + 2A* tatic A s t a tic) , Kt a tic(r)dV Wo cos(<? • s - c^)A* tatic (s) ,(77) 



9VN 
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FIG. 3: Bragg excitation of a Bogoliubov mode of wavevector ko and frequency uo for a finite temperature Bose condensed gas 
in a cubic box. The vertical dashed line at time t — it /ojo indicates the time after which the perturbation W is discontinued. 
Solid lines: evolution of the field amplitudes of the Bogoliubov modes with wavevectors ko = (12n/L,Q, 0) (upper curve) and 
— fco (lower curve) in the Wigner simulation after averaging over 100 realizations. Only the mode ko is excited resonantly 
by Bragg scattering. After the coherent excitation Bragg phase, the amplitudes of the two modes are damped. Dashed line: 
perturbative approach of subsection IIVB 21 The truncated Wigner approach and the perturbation theory give comparable 
results. N = 5 x 10 4 , k B T = 3/i, Twjo = 2.2/x, Wo = 0.175/x, fj, = 500h 2 /mL 2 . In the Wi gner simulation a grid with 22 points 
per dimension is used, so that J\f = 22 3 = 10648 <C N. In the perturbative approach a grid of 48 points per dimension is used 
to avoid truncation effects. The initial mean number of noncondensed particles is N — (No) ~ 5000. 



where Wo is non zero only during the excitation phase. In this equation ho = p 2 /2m + Wo cos(q ■ r — tot) is the 
one-body part of the Hamiltonian including the kinetic energy and the Bragg excitation potential, and Q projects 
orthogonally to the static condensate mode </>. The term of zeroth order in A sta tic is a source term which causes A sta tic 
to acquire a nonzero mean value during the evolution. The terms of first order in A stat ; c in (|77p describe the evolution 
in the static Bogoliubov approximation. Terms of second order provide the damping we are looking for. We project 
equation ((77|) over the static Bogoliubov modes fTj]) by using: 

Astatic (V) = ^ b k u k (r) + b* k v* k (r) (78) 

k^0 

with the mode functions u k (r) and v k (r) defined in (|74p . Terms nonlinear in A sta tic in (|77|) then correspond to an 
interaction between the Bogoliubov modes. 

We assume that the excitation phase is much shorter than the damping time of the coherently excited mode. As 
a consequence we can neglect in this phase the processes involving interaction among the Bogoliubov modes. Also 
in the action of the perturbation W we keep only the term acting on the condensate mode, that is the first term 
on the right hand side of ([77]) . which is yj (No) larger than the terms acting on the noncondensed modes. For the 
choice of parameters (1701) only the two modes with wavevectors k and — k are excited from the condensate by the 
perturbation W; the amplitudes of the field in these modes evolve according to 

ihjbo = huobo + ^A^^y(U + V )e- iuot (79) 

ih^b- = ^ 6_ + \/(]Vo)^([/o + Vo)e^ t . (80) 

By integrating these equations we realise that the mean amplitude (bo) grows linearly in time, since the mode is 
excited resonantly, while the mean amplitude (b—o) oscillates and vanishes at t = tt/ljo- 
After the excitation phase we include the second order terms that provide damping: 

ihjbo - <.,h, : + J2 A li b i b i + i A U + A o,i)Kt>j + J2(B id , + 11,.,,, + B h o,j)b*b* (81) 

with 

A U = ^^[U l (U J +V j )U k + (U l + V l )V ] U k + V 3 (U k + V k )V l ]5 l . ]+k (82) 
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/•',;./. = Qj^-ViiUj+VjW-u+t. (83) 

and where i,j,k denote momenta. The last terms with the £?'s in (fBTj) do not conserve the Bogoliubov energy and 
we can neglect them here for the calculation of the damping rate since we are going to use second order perturbation 
theory; we would have to keep them in order to calculate frequency shifts. In the terms with the A's we recognise two 
contributions: the term with A® . describes a Beliaev process where the excited mode can decay into two different 

modes while the term with A\ + A J , describes a Landau process where the excited mode by interacting with another 

mode is scattered into a third mode [29(. We introduce the coefficients b in the interaction picture 

bj = bj e ie i t/h (84) 



where tj is the Bogoliubov eigenenergy of the mode with wavevector j, and we solve ([81]) to second order of time- 
dependent perturbation theory to obtain: 

(bo(t) - 6o(0)) =s ~ £ A lMh + A li) 7 *( e ° ~ e * ~ ^X 1 + * + %)( fo o(0)) 

i,j 

~J2 EKo + A iif Itfa + d - - nj)(b o (0)) 

i,j 

-^2{Al+ a f I t (e + e - e Q+Q )(b*(0)b o (0)b o (0)) (85) 
where + represents the mode of wavevector 2ko and where 

It(v) = ! dre ivr ' n f T {v) (86) 
Jo 

f r {v) = [ d8e- w6 / h . (87) 
Jo 

The fij 's are the occupation numbers of the Bogoliubov modes in thermal equilibrium given by the Bose formula 

' (88) 



J e Cj/k B T _ i 

where €j is the energy of the Bogoliubov mode. In the language of nonlinear optics the last line in ([55)) describes a 
Xi effect or a second harmonic generation which can be important if the conservation of energy condition €2k a = 2e/c 
is satisfied and if the initial amplitude (&o(Q)) = y8 is large since one has 

(6*(0)6o(0)6o(Q)) = \(i\ 2 P + n 2/3 . (89) 

We have checked that the X2 effect is negligible for the low amplitude coherent excitations considered in the numerical 
examples of this paper: eo is larger than fi so that fco is n °t in the linear part of the Bogoliubov spectrum and therefore 
the second harmonic generation process is not resonant. By using the fact that: 

Re I t {v) = l -\f t { V )\ 2 = ^ sin 2 ^ = nMS t (v) (90) 

where 5t{v) converges to a Dirac delta distribution in the large t limit, we calculate the evolution of the modulus of 
the Bogoliubov mode amplitude 

|<6o(*)>|-K6o(0))| nt 



Y, A U A h + A li) S ^o -ei- e,)(l + m + fij) 



\(bo(0))\ 

~\ Yj A ^ + A U 2 S ^ + e^rii - n ) . (91) 
i,j 

This formula can be applied to a finite size box as it contains finite width <5's. By plotting equation (]91[) as a function 
of time we can identify a time interval over which it is approximately linear in time, and we determine the slope 
— 7perturb with a linear fit [3(|. Heuristically we then compare exp(— 7perturt>£) to the result of the truncated Wigner 
simulation, see figure [3] and we obtain a good agreement for this particular example [3l| . 

In the thermodynamic limit, when the Bogoliubov spectrum becomes continuous, the discrete sums in (1911) can be 
replaced by integrals and the finite width St is replaced by a Dirac 5 distribution. In this case an analytical expression 
for the damping rate can be worked out and we recover exactly the expression for the Beliaev and Landau damping 
rate obtained in the quantum field theory [H, [H, [34| . 
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3. Validity condition of the truncated Wigner approach 



We now investigate numerically the influence of the grid size on the predictions of the truncated Wigner simulation. 
The line with squares in figure 2] shows the damping rate obtained from the Wigner simulation, defined as the inverse 
of the 1/e half- width of \(b n (t))\, as a function of the inverse grid size 1/N. For small grids the results of the 
simulations reach a plateau close to the perturbative prediction 7 por turb- For large grids the damping rate in the 
simulation becomes significantly larger than 7 per turb- Since the perturbative prediction reproduces the known result 
for Beliaev-Landau damping, we conclude that the results of the truncated Wigner simulation become incorrect for 
large grid sizes. The reason of such a spurious damping appearing in the Wigner simulation for large N will become 
clear below. 
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FIG. 4: Damping rate of the coherent excitation in the Bogoliubov mode of wavevector fco = (12-w/L, 0, 0) and of frequency 
uo as a function of the inverse number of modes in the grid 1/Af for the Glauber-P and the Wigner distributions. Each disk 
represents the average over 100 realisations of the simulation and the lines are a guide to the eye. N = 10 5 , fcsT = 3fi, 
H — 500h 2 /mL 2 , so that Tiujq = 2.2/1, 7p ertur j, = 0.061mL 2 /?i, Wo = 0.0874/x. The damping rate is expressed in units of 
7perturb- Arrows indicate some values of e ma x/fcsT where e ma x is the maximal Bogoliubov energy on the grid. 

It is tempting to conclude from the perturbative calculation of subsection IIVB 21 that the validity condition of the 
truncated Wigner approach is dictated only by the condition N -C N. To check this statement we have performed a 
second set of simulations (not shown) for a particle number N reduced by a factor of two keeping the size of the box 
L, the chemical potential [i = Ng/L 3 and the temperature fixed. If the condition of validity of the truncated Wigner 
approach involves only the ratio N/N the plateaux in the damping time should start at the same value of N/N for 
the two sets of simulations. However this is not the case, and we have checked that on the contrary, the two curves 
seem to depend on the number of modes only. 

Another way to put it is that the condition to have agreement between the truncated Wigner simulation and the 
perturbation theory of section IIVB 21 is not (or not only) that the number of particles should be larger than the 
number of modes. There is in fact another "hidden" condition in the perturbative calculation which is the hypothesis 
that the occupation numbers of the Bogoliubov modes are constant during the evolution. In reality, even in absence 
of the Bragg perturbation, our initial state which reproduces the correct thermal distribution for the quantum Bose 
gas, is not stationary for the classical field evolution ([T]). The perturbative expression (f^T|) holds indeed in the limit 
N/N 3> 1, but the occupation numbers of the Bogoliubov modes, initially equal to the Bose formula fij, change in 
the course of the time evolution in the simulation and this affects the damping rate. This effect is neglected in the 
perturbative formula (|91[) and it is found numerically to take place on a time interval comparable to the damping 
time of the Bogoliubov coherent excitation as we show in figure [5] 

What it is expected to happen in the absence of external perturbation is that the classical field equation ([I]), in the 
three-dimensional cubic box geometry considered here, displays an ergodic behaviour leading to thermalisation of the 
classical field i/j towards its equilibrium distribution [111 . |l2j . In the regime where the noncondensed fraction is small 
and the number of modes is smaller than N, we can approximately view the classical field as a sum of Bogoliubov 
oscillators bk weakly coupled by terms leading to the nonlincarities in (|F2p . In the equilibrium state for the classical 
field dynamics we then expect the occupation numbers of the Bogoliubov modes to be given by the equipartition 
formula: 
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FIG. 5: Evolution of the squared amplitudes (blbk) of the classical field Bogoliubov modes multiplied by the corresponding 
Bogoliubov energy efe in the truncated Wigner simulation in the absence of the Bragg perturbation. We have collected the 
Bogoliubov modes in energy channels of width 2/x, so that the plotted quantity is the average among each channel of ek{blbk), 
with increasing energy from top to bottom at initial time t = 0. The thick horizontal line is the expected temperature T c i ass of 
the equilibrium classical field distribution as given by (|94[) . Parameters are: N = 5 ■ 10 4 , ksT = 3/i, fi = 500ft 2 /mL 2 and the 
vertical axis of the figure is in units of Ti 2 /mL 2 , where L is the cubic box size. The number of modes is 22 per spatial dimension 
so that the maximum Bogoliubov energy allowed on the grid is e ma x = 15.3/x. The averaging in the simulation is performed 
over 500 realisations. 

attributing a mean energy of kBT c \ ass to each of the Bogoliubov mode. The classical field equilibrium temperature 
Tciass can then be deduced from the approximate conservation of the Bogoliubov energy [35| : 



fc B T class = _L r ^ ejk (6*6*)(t = 0) 

fe/0 

1 ^ 



(93) 



exp(/?e fc ) - 1 2 

L_ V - (941 

- 1 ^ 2tanh(/W2l ' 1 ; 



TV - 1 ^2tanh(/3e fe /2) 

The thermalisation of the Bogoliubov modes to the new temperature T c i ass is nicely demonstrated in figure [5] One 
sees that tk{b* k bk) indeed converges to a constant value almost independent of k. From the fact that tanhx < x for any 
x > we deduce that the classical equilibrium temperature T c i ass is always larger than the real physical temperature 
T of the gas. In the regime fcgT ^> [i this 'heating' increases the squared amplitudes (b* k bk) of the modes of energy 
~ /i by a factor ~ T c \ aaa /T. Since the Landau damping rate is approximately proportional to the populations of these 
modes [32J, |33j, l34| the damping rate is increased roughly by a factor T c i ass /T, an artifact of the truncated Wigner 
approximation. 

It is clear that T c i ass will remain very close to T as long as the maximum Bogoliubov energy allowed in the simulation 
remains smaller than fc^T. One can indeed in this case expand (|94[) in powers of One has to expand the hyperbolic 
tangent up to cubic order to get a nonzero correction: 

Tciass ^ ^ , 1 (P e k) 



The absence of terms of order (3ek in (l95|) is a fortunate consequence of the noise added to the field in the Wigner 
representation. This added noise shifts the average (b* k bk){t = 0) by 1/2 with respect to the Bose formula. 

When the maximum Bogoliubov energy becomes much larger than fcgT we expect T c i ass to become significantly 
larger than T . This is illustrated in figure [6] obtained by a numerical calculation of the sum in (|94p for increasing 
grid sizes. We have also plotted in this figure the value that one would obtain for T c i ass in the absence of the added 
Wigner noise (i.e. in a Glauber-P approach), that is by removing the terms e^/2 in (f9"3"| . The Glauber-P distribution 
for the field ip in the sense of [3|| is given by 

v.- = aw + J2 bkUk + b >k ( 96 ) 
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where the bf, are chosen from a Gaussian distribution such that {b* k bk) — l/( ex P(/3 £ fc) — 1) an d the value of No is 
dictated by the normalisation condition H^H 2 = N. In this case T c i ass is always smaller than T, and deviates from T 
for smaller grid sizes, since the fortunate cancellation of the order f3ek obtained in ([95| does not occur anymore. We 
expect in this case a spurious reduction of the damping rate. We have checked it by evolving an ensemble of fields of 
the form ((96 j) with the Gross-Pitaevskii equation and we found that the damping rate is always smaller than half of 
the correct result even for the smallest grids that we tested, see the line with diamonds in figure |H 

2 | ■ , ■ , ■ , ■ , ■ — -, 
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FIG. 6: Equilibrium temperature Tciass of the classical gas as function of the maximum energy e max of the Bogoliubov modes 
on the momentum grid with the assumption of equipartition of the energy in the Bogoliubov modes. Circles: the initial field 
distribution is the Wigner distribution for the quantum gas at temperature T. Crosses: Glauber-P distribution defined in [36j], 
amounting to the removal of the added Wigner noise from the initial field distribution. The dashed lines are a guide to the eye. 
The number of momentum components along each dimension of space goes from 2 to 30 in steps of 2. The chemical potential 
is fi = 5007i 2 1 ml? and the temperature is ksT — 3/i. 



V. CONCLUSION 



We have considered a possible way of implementing the truncated Wigner approximation to study the time evolution 
of trapped Bose-Einstein condensates perturbed from an initial finite temperature equilibrium state. First a set of 
random classical fields ip is generated to approximately sample the initial quantum thermal equilibrium state of the 
gas, in the Bogoliubov approximation assuming a weakly interacting and almost pure Bose-Einstein condensate. Then 
each field ip is evolved in the classical field approximation, that is according to the time-dependent Gross-Pitaevskii 
equation, with the crucial difference with respect to the more traditional use of the Gross-Pitaevskii equation that 
the field ip is now the whole matter field rather than the field in the mode of the condensate. 

The central part of this paper is the investigation of the validity conditions of this formulation of the truncated 
Wigner approximation. 

For short evolution times of the fields ip the dynamics of the noncondensed modes, i.e. the components of the 
field orthogonal to the condensate mode, is approximately linear; we can then use the time-dependent Bogoliubov 
approximation, both for the exact quantum problem and for the truncated Wigner approach. A necessary condition 
for the truncated Wigner approach to correctly reproduce the quantum results is then 

N^>N/2 (97) 

where J\f is the number of modes in the Wigner approach and N is the total number of particles in the gas. This 
condition can in general be satisfied in the degenerate and weakly interacting regime without introducing truncation 
effects due to a too small number of modes. 

For longer evolution times the nonlinear dynamics of the noncondensed modes comes into play. When the classical 
field dynamics generated by the Gross-Pitaevskii equation is ergodic, e.g. in the example of a three dimensional gas 
in a cubic box considered in this paper, the set of Wigner fields ip evolves from the initial distribution mimicking the 
thermal state of the quantum gas at temperature T to a classical field equilibrium distribution at temperature T c i ass . 
Since noise is added in the Wigner representation in all modes of the classical field to mimic quantum fluctuations it 
turns out that r c i ass is always larger than T. If T c i ass deviates too much from T the truncated Wigner approximation 
can give incorrect predictions. For example we have found that the Beliaev-Landau damping of a Bogoliubov mode in 
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the box, taking place with a time scale comparable to that of the 'thermalisation' of the classical field, is accelerated 
in a spurious way as the classical field 'warms up'. A validity condition for the truncated Wigner approach in this 
long time regime is therefore 

|T class - T| « T. (98) 

This condition sets a constraint on the maximum energy of the Bogoliubov modes e max in the Wigner simulation: 
£max should not exceed a few kgT. More precisely one can use the following inequality to estimate the error [37| : 

\T ciass -T\ 1 (e|) 1 f e max \ 2 

T < 12(fc B T) 2 < 12 \k B Tj [ ' 

where (e|) is the arithmetic mean of the squares of all the Bogoliubov energies in the Wigner simulation. 

The fact that the initial set of Wigner fields is nonstationary under the classical field evolution could be a problem: 
the time-dependence of the observables could be affected in an unphysical way during the thermalisation to a classical 
distribution of the ensemble. To avoid this, we could start directly from the thermal equilibrium classical distribution 
[ijlEU' restricting to the regime e max < kgT. 

A remarkable feature of the Wigner simulation is that T c i ass deviates from T at low values of e m ax only quadratically 
in e max /kgT. This very fortunate feature originates from the added noise in the Wigner representation. It explains 
why for e max as high as 3.5 kgT the truncated Wigner approach can still give very good results for the Beliaev- 
Landau damping time (see Fig. [4|. In contrast, if we remove the Wigner added noise, in the so-called Glauber-P 
representation, or if we add more noise, in the so-called Q representation, T c i ass deviates from T linearly in e max /fcsT. 
In this case we expect that the condition of validity of the classical Gross-Pitaevskii equation will be that all modes 
in the problem must be highly occupied, resulting in the stringent condition e max < ksT. We therefore conclude 
that the Wigner representation is the most favorable representation of the quantum density operator with which to 
perform the classical field approximation. This fact, known in quantum optics for few mode systems, was not obvious 
for the highly multimode systems that are the finite temperature Bose gases. 

Still, condition (|98|) is a serious limitation of the truncated Wigner method for simula ting general ergodic three 
dimensional systems. One possibility to overcome this limitation is to proceed as in [38l |39| i.e. to treat the high 
energy modes as a reservoir, which leads to the inclusion of a stochastic term in the Gross-Pitaevskii equation. The 
advantage of this treatment is that the additional term has dissipative effects and thermalises the system to the 
correct quantum field thermal distribution in the stationary state as opposed to the classical one. However, one of the 
conceptual advantages of the truncated Wigner method and of classical field methods in general 0, [l(J, [111, [12] which 
we would like to keep is that apparent damping and irreversibility arise from the dynamics of a conservative equation 
(the Gross-Pitaevskii or nonlinear Schrodinger equation) as is the case in the original Hamiltonian equations for the 
quantum field. 

Laboratoire Kastlcr Brossel is a research unit of Ecole Normale Superieure and of Universite Pierre et Marie Curie, 
associated to CNRS. We acknowledge very useful discussions with Crispin Gardiner. This work was partially supported 
by National Computational Science Alliance under DMR 9900 16 N and used the NCSA SGI/CRAY Origin2000. 

APPENDIX A: BARE VS EFFECTIVE COUPLING CONSTANT 

In this appendix we describe how to adjust the potential V(r) defined on the grid in the simulation in order to 
reproduce correctly the low energy scattering properties of the true interatomic potential. 

We start with the Schrodinger equation for a scattering state (j>(r) of the discrete delta potential V(r) = (go/dV)S r o 
on the spatial grid of size L v and volume V: 

60(r) = ( P -<j) (r) + ^-Mr)8r,o (Al) 
\m J dv 

where m is twice the reduced mass and where (j>(0) is different from zero. We project this equation on plane waves of 
momentum k: 

~ m = ffL. *(°) , ( A2 ) 



where 4>(k) is the component of <f> on the plane wave e lk ' r /VV. Fourier transforming back gives 0(0); dividing the 
resulting equation by </>(0) leads to the quantization condition 

i = iy — f — . (as) 
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We define the effective coupling constant g e g in such a way that the energy of the lowest scattering state of the 
pseudopotential g e «S{r)d r (r •) in the box is the same as the energy of the lowest scattering state solution of (|A3|) . 

We now restrict ourselves to the case where the size of the box is much larger than the scattering length associated 
with g e fi . In this case the energy of the lowest scattering state for the continuous theory with the pseudopotential is 
very close to g e g/V, so that we can calculate g c g from the equation e = g e g/V. In this large box case, one can then 
check that the energy e is negligible as compared to h 2 k 2 /m except if k = 0. This gives 

1 "T V h 2 k 2 /m 

which allows us to adjust g in order to have g e fi — g = Airfi 2 a/ra where a is the scattering length of the true 
interatomic potential. 

The sum over k in the denominator can be estimated by replacing the sum by an integral over k and is found to 
be on the order of fc max ao where go = Anh ao/m and fc max is the maximum momentum on the grid, go is therefore 
very close to g e g when condition © is satisfied, so that we can set go ~ g e g = 9- In the opposite limit of a grid step 
size tending to zero one gets g e g — > 0, and we recover the known fact that a delta potential does not scatter in the 
continuous limit. We would have to increase go continuously up to infinity as the grid step size tended to zero, if we 
wanted to get a finite g e g in this limit. 




APPENDIX B: AN IMPROVED BROWNIAN MOTION SIMULATION 

A better choice for a and Y - In our previous work Q the drift matrix a and the noise matrix Y were the hyperbolic 
sine and cosine of C/(2kgT), which imposed a time step dt in the simulation which was exponentially small in the 
parameter e max /(fcsT), where e max is the largest eigenvalue of C allowed on the spatial grid of the simulation. We 
have now identified a choice that does not have this disadvantage: 

(Bl) 
(B2) 

where the projector Q is defined in (|25|) . With this new choice for a and Y both the friction matrix and the noise 
matrix are bounded from above by unity, which allows a much larger dt in the case e max > fc^T. To calculate the 
action of matrix a on the vector (tp±, ipj_) we write the hyperbolic tangent as: 

tanha; . . . 

ta,nhx = x = xF(x 2 ). (B3) 

x 

The function F{u) is then expanded on Chebyshev polynomials in the interval u £ [0, (e m ax/(2fcsT)) 2 ] and approxi- 
mated by a polynomial of a given degree, typically 15 for e max /(2A;BT) = 3 and 25 for e max /(2fesT) = 6, obtained by 
truncating a Chebyshev expansion of degree 50 |40l |. 

An improved integration scheme - Initially we set ip± = 0. Since the noise g?£ is Gaussian, and because the stochastic 
differential equation (|28p is linear, the probability distribution of ip± is guaranteed to be Gaussian at any step of the 
integration so that the issue of the convergence of the distribution to the correct steady state distribution ((2~l|) can 
be discussed in terms of the convergence of the covariance matrix of the distribution to its right steady state value. 
Two issues in particular should be addressed: the error introduced by the discretisation in time (finite time step dt 
of integration), and the error introduced by the integration over a finite time interval (approach to the steady state 
distribution) . 

We now explain how to face the first problem with an efficient integration scheme yielding an error on the steady 
state covariance matrix of the distribution scaling as dt 2 , rather than dt for the simple Euler scheme. In the numerical 
scheme the vector X = (ip±,ip*[) that stores the values of the field ^>± and of its complex conjugate on the discrete 
grid obeys the recursion relation: 

X[t=(n+i)dt\ = (1 — a num dt)X[ t=ndt ] + Y num I ^~ ndt ^ ) (B4) 



dt] 



with the initial condition Xu = o\ = 0. In this recursion relation the friction matrix a num and the noise matrix l^ um 
may differ from a and Y of the continuous stochastic differential equation (|28[) by terms linear in dt that remain to 
be determined in order to achieve an error scaling as dt 2 . 
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As we have already mentioned X\ t=n m is a Gaussian vector for any step n of the iteration so that its probability 
distribution is characterised by the covariance matrix = (XiXj), with indices i,j ranging from 1 to 2Af. From 
(|B4[) the covariance matrices are shown to obey the recursion relation: 

= (1 - a num dt)C^ (1 - al um dt) + ^r num F n t um . (B5) 

For a small enough time step dt this matrix sequence converges to a finite covariance matrix solving 

C (oo) = (1 - a num dt)C^(l - al um dt) + ^y num F n t um . (B6) 

We now try to choose the friction matrix and the noise matrix in order to minimise the deviation of C^°°^ from the 
desired value, which is the covariance matrix of the exact distribution (|2ip . equal to (2M dV)~ x . We look for a num 
and Y num differing from the theoretical values (|BllB2p by terms linear in dt, and leading to a covariance matrix 
different from the theoretical one by terms quadratic in dt: 

a num = 2M + a x dt (B7) 
Y num = pJ+Yidt (B8) 

CM = 5^ + (B9) 



Equation (|B6[) is satisfied up to order dt irrespectively of the choice of u\ , Y\. Requiring that equation (|B6[) is satisfied 
up to order dt 2 leads to the condition 

- ai ±--±-a 1+Yl ( Q °) + ( Q ° W + M = 0. (BIO) 
4M 4M \0Q*J \0QJ 

A particular solution of this equation is provided by a% — and Y\ = = —M/2. Our improved integration scheme 
is therefore 

a num = 2M (Bll) 

Y num = ( Q q M- l -Mdt. (B12) 

The analysis of the recursion relation (|B5|) is easily performed for our improved integration scheme (|B1 1IF312|) since 
Q^num, Q^numi ^num and hence C^ n ' are polynomials of M and commute with M. As a consequence C*- 00 ^ also commutes 
with M. 

Let us first estimate the deviation of C^ 00 -* from the exact covariance matrix (2M dV) -1 : 



£f(oo) _ 



1 - (1 - a num di) J (B13) 



1 



2MdV 



l + ^M 2 + 0(dt 3 ) 
4 



(B14) 



Because M is bounded from above by unity we take in practice dt = 1/8 so that the error is less than 0.5 percent. 
Let us finally estimate the convergence time of the covariance matrices. The recursion relation (|B5[) can be rewritten 

as 

C (n+1) - C (oo) = (1 - a mm dtf [C (n) - C (oo) j (B15) 

so that the relative deviation of C^ n ' from its asymptotic value evolves as (1 — 2M min dt) 2n where M m i n is the smallest 
eigenvalue of M, that can be evaluated along the lines of [J]. We choose the number of time steps n so that the 
relative deviation of from C^ 00 -* is less than 0.5 percent. 
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APPENDIX C: MOMENTS OF N OF A HARMONICALLY TRAPPED IDEAL BOSE CONDENSED GAS 



We explain how to calculate the approximate expressions (|42[) for the moments of the number of condensed particles 
for an ideal Bose gas in an isotropic harmonic potential of frequency u> in the temperature regime kgT 3> tuv and 
in the Bogoliubov approximation. The calculation of the moments involves sums over the excited harmonic levels, 
see (|4Tj) . By using the known degeneracy of the harmonic eigenstate manifold of energy nhio above the ground state 
energy the calculation reduces to the evaluation of sums of the form 



■We) = Y. 



^ (exp(ne) - 1)9 



(CI) 



where e = Hw/ksT is tending to zero, and the exponents p and q are positive integers. 

First case: q — p > 1: In the limit e — > the sum is dominated by the contribution of small values of n. Replacing 
exp(en) — 1 by its first order expression we obtain: 



1 °° 1 1 

n—l 



(C2) 



where C( a ) = S n >i l/ n< * i s the Riemann Zeta function. 

Second case: q — p < 1: In the limit e — > the contribution to the sum is dominated by large values of n. We then 
replace the discrete sum by an integral over n from 1 to +oo. Taking as integration variable u = en we arrive at 



S p , q (e) 



1 



+oo 



du ■ 



■IP 



We can take the limit e 



e p+1 h (exp(tt) - 1)9 ' 

in the lower bound of the integral since q — p < 1: 

1 



In 



(C3) 



(C4) 



To calculate the resulting integral I Ptq we expand the integrand in series of exp(— u) and integrate term by term over 
u: 



du 



P 



(exp(w) - l) 9 



^(k + qf+i k\(q-l)\ 



(C5) 



which can be expressed in terms of the Riemann Zeta function, e.g. ^2,2 = 2(C(2) — C(3))- 

Third case: q — p = 1 : In the limit e both the small values of n and the large values of n contribute to the sum. 
We introduce a small parameter i/ < 1 that will be put to zero at the end of the calculation. For the summation 
indices n < v je we keep a discrete sum and we approximate each term of the sum by its first order expression in 
e, which is correct as ne < v <C 1. For the summation indices n > vje we replace the sum by an integral, which 
is correct in the limit e — > for a fixed v, since we then recognise a Riemann sum of a function with a converging 
integral. This leads to 



1 



u/e 



du 



In the limit e — > the discrete sum is approximated by 



£i^log(*Ve) 



(exp(u) - 1)p +1 



+ 7 



(C6) 



(C7) 



where 7 is Euler's constant. In the integral we remove and add l/(exp(u) — 1) to the integrand in order to get a 
convergent integrand which facilitates the calculation of the v — > limit. The integral of l/(exp(u) — 1) can be 
calculated explicitly from the primitive log(l — exp(— u)) so that in the small v limit 



+00 



du 



(exp(u) - 1) 



1 



cxp(— v) 
log 1/ + J p 



1 



(exp(u) — 1)p +1 exp(u) — 1 



(C8) 
(C9) 
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where 



+00 



du 



(exp(w) — exp(w) — 1 



(CIO) 



The — log v term coming from the integral compensates the log v term coming from the sum in (|C7[) so that in the 
limit v — * we get the ^-independent estimate 



1 



Sp,p+1 ^ -f^ [~ log e + 7 + J"p] • 



(Cll) 



The quantity J p for p > can be calculated from a recursion relation obtained in the following way: we use the 
identity 



exp(uj 



(exp(w) - (exp(w) - 1)p (exp(w) - 



(C12) 



The first term of the above expression leads to an integral already calculated in (|C5|) and called I PtP . We then integrate 
the second term of the above expression by parts, taking the derivate of u p with respect to u. This finally leads to 



Jp Jp—i T Ip,p- 



(C13) 



We get in particular J x = 1 - C(2) and J 2 = 3/2 - 3<C(2) + 2((3). 

Finally we collect the approximations for the S p ^ q relevant for the calculation of the skewness of the number of 
condensed particles No in ID, 2D, 3D: 



So, 



Si, 



Tog e + 7 



C(2) 



So, 2 

<Sl,2 



C(2) 



So,: 



C(3) 



log(e)+7 + l-C(2) q C(2) 



(C14) 



S2.1 



2C(3) 
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2C(2) - 2C(3) 



S-. 



2,3 



log e + 7 + J 2 



APPENDIX D: EQUATIONS OF THE NUMBER CONSERVING BOGOLIUBOV APPROACH 



In this appendix we give the equations of motion for the operator A and for (j)^ 1 (r) from [B| . The evolution equation 
for A is: 



ihd t 



A(r,t) 



with C given by ([24]) . The evolution equation for (p]_'(r) is: 

d 



= £(t) 



ih— - C(t) 
dt K ' 



A(r,t) 
At (r,t) 



Q(t)S(t) \ 
Q*(t)S*(t) ) 



(Dl) 



(D2) 



where 



S(r) = -giV|0(r)| 2 0(r)(l+^^At( s )A( s )) 

S 

+ 2gN4>(r)(A 1i (r)A(r)} + gN<f>*(r){A(r)A(r)) 

- gN^dVl^s) | 2 <[At( s )0( s ) + A( S )0*( S )j A(r)> 



(D3) 
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APPENDIX E: TRUNCATED WIGNER APPROACH IN THE BOGOLIUBOV REGIME 

In this appendix we demonstrate the equivalences (|60H63|) . For convenience we change in this appendix the phase 
reference of the field ip which now evolves according to 

ihdttp = [p 2 /2m + U(r, t) + g\iP\ 2 -/j]tp (El) 

where /j, is the chemical potential in the time-independent Gross-Pitaevskii equation for the condensate wavefunction 

m- 

1. Identification of the pure condensate wavefunction 

At t = equation (|60|) is satisfied. By keeping only terms of order y/~N in (|E1|) , in the limit (j4~4"|) , we obtain 

ihd t ^ 0) = (h Q + g\ V> (0) | 2 - a0^ (o) (E2) 
where ho is the one-body part of the Hamiltonian. This shows that (f60|) holds at all times. 

2. "Orthogonal-orthogonal" contribution 

We wish to prove (|62p . To this aim we expand A and tp± over the Bogoliubov modes: 

A = ]T ku k + b\v* k (E3) 

k 

^ =J2hu k + b* k v* k (E4) 

k 

At t = the same mode functions u k and v* k appear in the expansions of A and tp± . We wish to show that 

(|E3IIE4j) hold at any time, or equivalently that A and have the same equations of motion. If we keep only 
terms of order 0(1) in (|E1[) we get 

m {$l)= C °r(^») < E5) 

where Cqp is the usual Bogoliubov operator obtained from (j2"4")) by eliminating all the projectors. By using the 
fact that 



and 



with the matrices 



ViM _ { Q o ) ^« 

i I Q* { ^> 



£W<t> \ = (V \ 



(E6) 



(E7) 



V r , s = dV <t>{r)ct>* (s) Q riS = S r , s - dV^(r)0* (s) (E8) 



we get 



ihd t 



V ^ ^ ) + (e + ^ } I o s -J ( -,*|*| V J (E9) 

= dy^gAr|^(r)| 2 [0*(r)^W(r)+V (1) *(r)^(r)]. (E10) 

r 

The fact that the derivative of £W is purely imaginary and the initial condition £W = guarantee that 
+ = for all times, which proves that A and have the same equations of motion. At all times 
we then have 

(At(s)A(r)) = Uk(r)u* k (s)(blb k ) + v* k (r)v k (s)(b k bl) (Ell) 
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and 

(VA^OO^M) = (At(s)A(r)) + ^ Uk (r)ut( S )-v* k (r)v k ( S ) (E12) 

where the amplitudes b k are time-independent and the u k , v k are time-dependent modes evolving according to 

ihth ( "* ) = £ ( "* ) ■ (E13) 

By using the decomposition of unity, equation (61) of reference 

u k(r)<(s) ~ v* k (r)v k (s) = ^Qr, s (E14) 

k 

we prove (|6"2")l . 

3. "Parallel-parallel" contribution 

We wish to prove (j6Tj) . We use the fact that (dVJ2 r W 7 ")! 2 ) ^ s a constant of motion order by order in 1/yN. 

^-N = (E15) 



To order V N we get 



To order iV° we get 



df 

which we verified directly in (|E10p . To order 1 / y/N we get 



d =0 (E16) 



(e (2) + e (2) *> + (i(e (1) i 2 > + (^EiA 1) M^ 



(E17) 



Using ([Hll) we then obtain 



(? (2) + £ (2) *> + (|(e (1) | 2 > + (^> + ^tt^ = constant . (E18) 



At t = from (ISSi (l56l) we deduce 



A/"- 1 , 

constant = — - — (E19) 

so that at any time 

<C (2) +£ (2) *) + <|(C (1) | 2 ) = ~(SN) . (E20) 



Note that without the approximation in [3] we would have at t = constant = 4^ and as a consequence 
(£(2) +£(2)*) + (|(£ (1) | 2 ) = ~{5N) + \. The contribution of the 1/2 compensates exactly the term -\4>*{s)(j){r) 
in (|55|) . We neglect here this contribution. 



4. Term "parallel-orthogonal" 



The last step consists in proving (p)3"|) . We first remark that at t = ("01 ) = ^, and for linearity reasons 

(^l ) = at all times. At i = (|1)3")) is satisfied by construction. We then have to deduce the equation of 
motion for 

( X ) = (£«*^W+^(2)) (E21) 



2G 



and show that it coincides with the equation of motion for (j)^ . By keeping only terms of order 1/yfN in (|E1[) 
we get 

iha ( ^ (2) \-r ( ^ ^ + ( 9N[Fi> {l)2 + 20|V (1) | 2 ] \ (F22) 

%Wt { ^> )~^[ ^> ) + { -gN[MM« + 2r|V (1) | 2 ] ) ' ( } 

With a calculation analogous to the one we performed to obtain the derivative of (tf)± ,ipj^*), using (|E18[) to 
eliminate and replacing ipW by £^(j> + ip± \ we obtain: 

iUO t ( C )=c( C ) - m { ^t^t ) (^23) 



-gNQ*[2\4, { l ) \ 2 cj ) + 2^cf 2 ^l ) * + </>*(^) 



(E24) 



In particular, we find that the terms involving |£^| 2 disappear because (£^) 2 = — |C | 2 - By using (|E9|) and 
()E10|) we can calculate the derivative of (%): 



with 



JZ(r) = ~(5N)gN\<f > (r)\ 2 <t>(r) + 2gN<t>(r)l(tiA)~±\<t ) (r)\ 2 ^ 



+ gNt^A 2 )-gNl^<jy(rmr)\ 2 +dVj2 W*)| 2 <[ At («M*) + ^*(s)A( a )]A(r)> | (E26) 

which is identical to (ID3|) . except for the contribution of the term 1/2 neglected in Q as discussed after (|E20|) . 
In order to obtain (|E26[) we used the identity (f6"2"]) and the fact that all terms proportional to 0(r) are killed by 
the projector Q in (|E25|) . Summarising, (|E25|) and (|E26|) together with (ip£ ') = prove (I63)) . 



APPENDIX F: EQUATION FOR THE NONCONDENSED FIELD IN THE WIGNER APPROACH 



In the truncated Wigner approach, we define the field A cx (r) = a^j_(r)/ V N where </> is at this stage an arbitrary 
wave function normalised to unity, tp± is the component of tp orthogonal to <p, and a^, is the coefficient of ip along <f>. 
When ip solves the time-dependent Gross-Pitaevskii equation, the equation of motion for A cx is given by: 



ih 



dA c > 
dt 



Rk{r, s) 



s k=0 



N dt 



where we have collected the terms of the same power in A e 
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R 3 (r,s) = gN 



N 



i? 4 (r, s) 



2gN^-\^s)\ 2 A* cx (s)A cx (s)A cx (r) 
-gN ( A£(s)A ex (s)A ex (r)4,(s) 



N 



(F2) 



where = a^a^, h — p 2 /2m + U(r,t) is the one-body part of the Hamiltonian and Q r . s = 5 rjS — dV<p(r)<p* (s) 

projects orthogonally to <fi. In the case of a uniform wavefunction <f>(r) = 1/L 3 / 2 we have the following simplifications: 
(i) dt<f> is equal to zero, (ii) the constant terms like |0(s)| 2 0(s) are killed by the projectors, (hi) for terms having a 
vanishing spatial sum, can be replaced by S r ^ s , (iv) the sum over s of t/j±( s ) an d therefore of A ox (s) is zero. For 
this value of tfi, A cx coincides with A sta tic defined in (172)) and N$ is equal to iVo of equation (|77p . 
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